#Copyright (C) Nial Peters 2012
#
#This file is part of AvoScan.
#
#AvoScan is free software: you can redistribute it and/or modify
#it under the terms of the GNU General Public License as published by
#the Free Software Foundation, either version 3 of the License, or
#(at your option) any later version.
#
#AvoScan is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#GNU General Public License for more details.
#
#You should have received a copy of the GNU General Public License
#along with AvoScan.  If not, see <http://www.gnu.org/licenses/>.

import numpy
import datetime
from avoscan import processing
import sys
import scipy
from scipy.interpolate import interp1d
import math
import os.path
import pylab as plt
from ica import cPlume
from avoscan.processing import date2secs, fit_gaussian, calc_plume_height_from_geometry

def test_calc_plume_height():
    
    f1 = '/home/njp39/volcanology/CHILE/Chile_2012/Lascar/25Feb2012/DOAS/Scanner/2ndposscan1/Processed/20120225/25feb_2ndpossscan1__wangle.txt'
    f2 = '/home/njp39/volcanology/CHILE/Chile_2012/Lascar/25Feb2012/DOAS/Scanner/scan1/Processed/20120225/25feb_scan1_wangle.txt'
    
    #find the two scans which are closest in time to each other
    scans1 = [s for s in fullscans(f1)]
    scans2 = [s for s in fullscans(f2)]
    
    times1 = [min(s['times']) for s in scans1]
    times2 = [min(s['times']) for s in scans2]
    print times1[:10]
    print times2[:10]
    
    
    all_times = [sorted(times1)[5]] + times2
    all_times.sort()
    print "\n\n\n"
    print all_times[:10]
    
    delta_t = []
    for i in range(len(all_times)-1):
        delta_t.append(abs(all_times[i]-all_times[i+1]))
    
    print "min start time difference = ",min(delta_t)
    
    min_t_idx = delta_t.index(min(delta_t))
    print "min t idx = ", min_t_idx
    
    t1 = all_times[min_t_idx]
    t2 = all_times[min_t_idx + 1]
    
    try:
        scan1_idx = times1.index(t1)
        scan2_idx = times2.index(t2)
    except ValueError:
        scan1_idx = times1.index(t2)
        scan2_idx = times2.index(t1)
    
    print "scan1 idx = ", scan1_idx
    print "scan2 idx = ", scan2_idx
    
    for j in range(len(scans1)):
        if min(scans1[j]['times']) == t1:
            i = scans1[j]
            break
        elif min(scans1[j]['times']) == t2:
            i = scans1[j]
            break
    

    
    s1 = ProcessedScan(i["angle"].copy(), i["times"].copy(),i["so2_col_amount"].copy(), i["angle_error"].copy(), i["so2_col_amount_error"].copy(), i["int_times"].copy(), i["coadds"].copy(), 0.0, 0.0, 0.0, 0.0, 0.0)
    for j in range(len(scans2)):
        if min(scans2[j]['times']) == t1:
            i = scans2[j]
            break
        elif min(scans2[j]['times']) == t2:
            i = scans2[j]
            break
    s2 = ProcessedScan(i["angle"].copy(), i["times"].copy(),i["so2_col_amount"].copy(), i["angle_error"].copy(), i["so2_col_amount_error"].copy(), i["int_times"].copy(), i["coadds"].copy(), 0.0, 0.0, 0.0, 0.0, 0.0)
    
    
    print "Scan 1 at time = ", min(s1.times)
    print "Scan 2 at time = ", min(s2.times)
    
    print calculate_plume_height(s1, s2, (626674.789, -2579874.177, 4135),(625854.215, -2579206.115,3968.0))
    

def test_calc_plume_height2():
    
    f1 = '/home/njp39/volcanology/CHILE/Chile_2012/Lascar/29Feb2012/DOAS/Scanner/scan1/Processed/20120229/29feb_scan1_wangle.txt'
    f2 = '/home/njp39/volcanology/CHILE/Chile_2012/Lascar/29Feb2012/DOAS/Scanner/junctionscan_2/Processed/20120229/29feb_junctionscan_2_wangle.txt'
    
    #find the two scans which are closest in time to each other
    scans1 = [s for s in fullscans(f1)]
    scans2 = [s for s in fullscans(f2)]
    
    times1 = [min(s['times']) for s in scans1]
    times2 = [min(s['times']) for s in scans2]
    print times1[:10]
    print times2[:10]
    
    
    all_times = [sorted(times1)[2]] + sorted(times2)
    all_times.sort()
    print "\n\n\n"
    print all_times[:10]
    
    delta_t = []
    for i in range(len(all_times)-1):
        delta_t.append(abs(all_times[i]-all_times[i+1]))
    
    print "min start time difference = ",min(delta_t)
    
    min_t_idx = delta_t.index(min(delta_t))
    print "min t idx = ", min_t_idx
    
    t1 = all_times[min_t_idx]
    t2 = all_times[min_t_idx + 1]
    
    try:
        scan1_idx = times1.index(t1)
        scan2_idx = times2.index(t2)
    except ValueError:
        scan1_idx = times1.index(t2)
        scan2_idx = times2.index(t1)
    
    print "scan1 idx = ", scan1_idx
    print "scan2 idx = ", scan2_idx
    
    for j in range(len(scans1)):
        if min(scans1[j]['times']) == t1:
            i = scans1[j]
            break
        elif min(scans1[j]['times']) == t2:
            i = scans1[j]
            break
    

    
    s1 = ProcessedScan(i["angle"].copy(), i["times"].copy(),i["so2_col_amount"].copy(), i["angle_error"].copy(), i["so2_col_amount_error"].copy(), i["int_times"].copy(), i["coadds"].copy(), 0.0, 0.0, 0.0, 0.0, 0.0)
    for j in range(len(scans2)):
        if min(scans2[j]['times']) == t1:
            i = scans2[j]
            break
        elif min(scans2[j]['times']) == t2:
            i = scans2[j]
            break
    s2 = ProcessedScan(i["angle"].copy(), i["times"].copy(),i["so2_col_amount"].copy(), i["angle_error"].copy(), i["so2_col_amount_error"].copy(), i["int_times"].copy(), i["coadds"].copy(), 0.0, 0.0, 0.0, 0.0, 0.0)
    
    
    print "Scan 1 at time = ", min(s1.times)
    print "Scan 2 at time = ", min(s2.times)
    
    print calculate_plume_height(s2,s1,(621326.96258108818, -2582703.8573431349,3860), (622178.26182869286, -2583487.0817963425, 3927))
    
    
    
    


def calculate_plume_height(scan1, scan2, pos1, pos2):
    """
    Calculates the horizontal and vertical distances of the plume from the two scan
    locations. 
    """
    #TODO - function under construction - remember to update the doc string
    if not isinstance(scan1,ProcessedScan) or not isinstance(scan2,ProcessedScan):
        print type(scan1)
        raise TypeError, "Scan arguments must be ProcessedScan objects"
    
    if type(pos1) != tuple or type(pos2) != tuple:
        raise TypeError, "Position arguments must be an (x,y,z) tuple"
    
    #fit gaussians to the two scans to determine the plume position
    scan1._ProcessedScan__calculate_bkgd()
    scan2._ProcessedScan__calculate_bkgd()
    plume1 = scan1.get_bkgd_subtracted_plume_points()
    plume2 = scan2.get_bkgd_subtracted_plume_points()
    
    plt.plot(scan1.angles, scan1.col_amount, "r-")
    plt.plot(scan2.angles, scan2.col_amount, "g-")
    plt.show()
    #plt.plot(plume1[1], plume1[2], "r-")
    #plt.plot(plume2[1], plume2[2],"g-")
    #plt.show()
    
    plume_peak1 = fit_gaussian(plume1[1], plume1[2]) #[1] = angles, [2] = col_amounts
    plume_peak2 = fit_gaussian(plume2[1], plume2[2]) #[1] = angles, [2] = col_amounts
    
    print "Angle1 = ",plume_peak1.mean
    print "Angle2 = ",plume_peak2.mean
    
    return calc_plume_height_from_geometry(pos1, pos2, plume_peak1.mean, plume_peak2.mean)



class ProcessedScan:
    def __init__(self, angles, times, col_amounts, angle_err, col_err, int_times, coadds, h_dist, h_dist_err, v_angle, v_angle_err, int_method, bkgd_fit_order=1, bkgd_threshold=3.0):
        self.__is_calculated = False
        
        self.angles = angles
        self.times = times    
        self.col_amount = col_amounts
        self.angle_err = angle_err
        self.col_err = col_err
        self.int_times = int_times
        self.coadds = coadds
        self.h_dist = h_dist
        self.h_dist_err = h_dist_err
        self.v_angle = v_angle
        self.v_angle_err = v_angle_err
        self.integration_method = int_method
        
        #background fitting parameters
        self.bkgd_fit_order = bkgd_fit_order
        self.bkgd_threshold = bkgd_threshold
        self.manual_bkgd = False
        
        #error calc parameters
        self.error_calc_iters = 1 #start by not doing any error analysis
    
    
    def angle2time(self, angle):
        if angle < numpy.min(self.angles) or angle > numpy.max(self.angles):
            return None
        
        c = numpy.rec.fromarrays([self.angles, self.times])
        c.sort() #sort based on angles

        interp = interp1d(c.f0, numpy.array([date2secs(d) for d in c.f1]))
       
        return datetime.datetime.fromtimestamp(interp(angle))
        
    
    def get_bkgd_subtracted_plume_points(self):
        """
        Returns background subtracted col_amounts, times, angles etc....just
        for points in the plume (i.e. between the background limits that have been set)
        """
        
        #TODO - what if this stuff hasn't been calculated yet?
        plume_mask = numpy.logical_not(self.bkgd_mask)
        plume_idxs = numpy.where(plume_mask)        
        bkgd_subtracted_ydata = self.col_amount - self.bkgd_func(self.angles)
        return self.times[plume_idxs], self.angles[plume_idxs], bkgd_subtracted_ydata[plume_idxs], self.angle_err[plume_idxs], self.col_err[plume_idxs], self.coadds[plume_idxs], self.int_times[plume_idxs]
        
                
    def calculate(self):
        #choose the calculation method based on the selected integration method
        print self.times[0]
        if self.integration_method == "Trapezoid Rule":
            return self.__calc_trapezoid_method()
        elif self.integration_method == "Cylinder Fit":
            return self.__calc_cylinder_fit_method()
        else:
            raise ValueError, "Unknown integration method \'"+str(self.integration_method)+"\'"
        
    
    def __calc_trapezoid_method(self):
        self.__calculate_bkgd()
        calculate_ica_trapz(self)
        return self
    
    def __calc_cylinder_fit_method(self):
        self.__calculate_bkgd()
        self.__estimate_plume_radii_extremes()
        calculate_ica_cyl(self, self.error_calc_iters)
        self.__is_calculated = True
        return self
    
    
    def plot_radius_fit(self, ax):
        plume_mask = numpy.logical_not(self.bkgd_mask)
        plume_idxs = numpy.where(plume_mask)
        if self.ica == 0:
            ax.clear()
            return
        bkgd_subtracted_ydata = self.col_amount[plume_idxs] - self.bkgd_func(self.angles)[plume_idxs]
        ax.plot(self.angles[plume_idxs]+self.angle_err[plume_idxs], bkgd_subtracted_ydata + self.col_err[plume_idxs], 'g-', linewidth=0.5)
        ax.plot(self.angles[plume_idxs]+self.angle_err[plume_idxs], bkgd_subtracted_ydata- self.col_err[plume_idxs], 'g-', linewidth=0.5)
        ax.plot(self.angles[plume_idxs]-self.angle_err[plume_idxs], bkgd_subtracted_ydata + self.col_err[plume_idxs], 'g-', linewidth=0.5)
        ax.plot(self.angles[plume_idxs]-self.angle_err[plume_idxs], bkgd_subtracted_ydata - self.col_err[plume_idxs], 'g-', linewidth=0.5)
        
        #plot the ellipse fit at higher angle resolution
        high_res_ang = numpy.linspace(numpy.min(self.angles[plume_idxs]-self.angle_err[plume_idxs]), numpy.max(self.angles[plume_idxs]+self.angle_err[plume_idxs]), 1000)
        
        for radius, theta0, scaling in self.plume_radii:
            ax.plot(high_res_ang,scaling*length_of_scan_path(radius, theta0, self.h_dist, self.v_angle, high_res_ang), 'r-')
    
    
    def plot_bkgd_fit(self, ax):
        
        ignored = zip(*self.bkgd_ignored_pts)
        if ignored:
            ax.errorbar(ignored[0], ignored[1], ignored[3], ignored[2], 'r+')
        
        #subplot_calc.errorbar(not_bkgd_angles,not_bkgd_col_amounts, not_bkgd_col_err, not_bkgd_angle_err,'g+')        
        ax.errorbar(self.angles[self.bkgd_idxs], self.col_amount[self.bkgd_idxs], self.col_err[self.bkgd_idxs], self.angle_err[self.bkgd_idxs], 'b+')
        ax.plot(self.angles, [self.bkgd_func(i) for i in self.angles], 'm-', linewidth=2)
   
   
      
    def __calculate_bkgd(self):
        
        if not self.manual_bkgd:
            #Do the Gaussian fitting in order to estimate the width of the plume and it's position
            if len(self.angles) != len(self.col_amount):
                print "lengths do not match"
            try:
                gfit_params = processing.fit_gaussian(self.angles, self.col_amount)
                #mask the angles and col_amounts array to only expose the out of plume sections
                self.bkgd_lowlim = gfit_params.mean - math.fabs(self.bkgd_threshold * gfit_params.sigma)
                self.bkgd_highlim = gfit_params.mean + math.fabs(self.bkgd_threshold * gfit_params.sigma)
            
                if self.bkgd_lowlim > self.bkgd_highlim or self.bkgd_lowlim < min(self.angles[0],self.angles[-1]): #remember scans can be in either direction
                    #this quite often happens if there are insufficient points in the plume to effectively
                    #fit a gaussian to it. As a hack - we just reset the lims to include all of the data
                    # there is probably a smarter way to do this
                    print "time = ", self.times[0], "mean = ", gfit_params.mean, "sigma = ", gfit_params.sigma
                    self.bkgd_lowlim = gfit_params.mean
                    self.bkgd_highlim = gfit_params.mean
                
            except RuntimeError:
                #failed to fit gaussian - say that whole spectrum is background
                self.bkgd_lowlim = 90
                self.bkgd_highlim = 90

        #don't allow the background to be bigger than the scan data
        if self.bkgd_lowlim < min(self.angles[0], self.angles[-1]):
            self.bkgd_lowlim = min(self.angles[0], self.angles[-1]) 
        
        if self.bkgd_highlim > max(self.angles[0], self.angles[-1]):
            self.bkgd_highlim = max(self.angles[0], self.angles[-1])        
        
        is_lower_bkgd = self.angles < self.bkgd_lowlim
        is_higher_bkgd = self.angles > self.bkgd_highlim
        self.bkgd_mask = numpy.logical_or(is_lower_bkgd, is_higher_bkgd)
        bkgd_idxs = numpy.where(self.bkgd_mask)
        
        #fit a polynomial through the background
        ignored_pts = []
        print "bkgd_order = ",self.bkgd_fit_order
        while True:
            bkgd_angles = self.angles[bkgd_idxs]
            bkgd_col_amounts = self.col_amount[bkgd_idxs]
            bkgd_angle_err = self.angle_err[bkgd_idxs]
            bkgd_col_err = self.col_err[bkgd_idxs]
            
            polyfit_params = numpy.polyfit(bkgd_angles, bkgd_col_amounts, self.bkgd_fit_order)
            bkgd_function = numpy.poly1d(polyfit_params)
                 
            #calculate the residuals - ignoring those that are within error
            resid = find_residuals(bkgd_angles, bkgd_col_amounts, bkgd_angle_err, bkgd_col_err, bkgd_function)
            
            if len(numpy.nonzero(resid)[0]) == 0:
                break
            
            #remove the point with the greatest residual from the fit
            idx_to_skip = bkgd_idxs[0][numpy.argmax(resid)]
            ignored_pts.append((self.angles[idx_to_skip], self.col_amount[idx_to_skip], self.angle_err[idx_to_skip], self.col_err[idx_to_skip]))
            bkgd_idxs = (bkgd_idxs[0][numpy.where(bkgd_idxs[0] != idx_to_skip)],)
        
        #store the calculated background for all angles
        self.bkgd_func = bkgd_function
        self.bkgd_ignored_pts = ignored_pts
        self.bkgd_idxs = bkgd_idxs

        
    def __estimate_plume_radii_extremes(self):
        
        plume_mask = numpy.logical_not(self.bkgd_mask)
        plume_idxs = numpy.where(plume_mask)
        
        bkgd_subtracted_ydata = self.col_amount - self.bkgd_func(self.angles)
        
        #do the calculation for the full range of errors
        tmp = [(self.angles[plume_idxs] + self.angle_err[plume_idxs], bkgd_subtracted_ydata[plume_idxs] + self.col_err[plume_idxs]),
               (self.angles[plume_idxs] + self.angle_err[plume_idxs], bkgd_subtracted_ydata[plume_idxs] - self.col_err[plume_idxs]),
               (self.angles[plume_idxs] - self.angle_err[plume_idxs], bkgd_subtracted_ydata[plume_idxs] + self.col_err[plume_idxs]),
               (self.angles[plume_idxs] - self.angle_err[plume_idxs], bkgd_subtracted_ydata[plume_idxs] - self.col_err[plume_idxs])]

        self.plume_radii = []
        for ang, so2 in tmp:
            if len(so2) < 5:
                self.ica = 0
                self.ica_err = 0
                min_time = min(self.times)
                max_time = max(self.times)
                self.ica_time = min_time + ((min_time - max_time) / 2)
                self.ica_time_err = abs(self.ica_time - min_time)
                return

            #first of all we estimate the plume radius
            # guess some fit parameters - we use the width at half height as an estimate of radius
            ymax = so2.max()
        
            angle_at_half_height = ang[numpy.abs(so2 - (ymax / 2.0)).argmin()]
               
            theta0_guess = numpy.average(ang,weights=numpy.abs(so2))
            
            variance = numpy.dot(numpy.abs(so2), (ang-theta0_guess)**2)/numpy.abs(so2).sum()  # Fast and numerically precise
            
            r_guess_angle = math.sqrt(variance)
            r_guess = self.h_dist * math.tan(math.radians(r_guess_angle))
    
            scaling_guess = ymax / (2.0 * (r_guess / math.cos(math.radians(self.v_angle))))
        
            #put guess params into an array ready for fitting
            p0 = numpy.array([r_guess, theta0_guess, scaling_guess])
        
            #define the path length function and associated error function
            fitfunc = lambda p, x: p[2] * cPlume.get_path_length(p[0], p[1], self.h_dist, self.v_angle, x)
            errfunc = lambda p, x, y: fitfunc(p, x) - y
           
            # do the fitting
            p1, success = scipy.optimize.leastsq(errfunc, p0, args=(ang, so2))
        
            if success not in (1, 2, 3, 4):
                raise RuntimeError, "Could not fit cylindrical plume to data."

        
            self.plume_radii.append((p1[0], p1[1], p1[2]))
            



def length_of_scan_path(r, theta0, h_dist, v_angle, theta):
    """
    By assuming a cylindrical plume, this function returns the path length through the
    plume for each scan angle. For a cylindrical plume, the path through the plume
    will always be an ellipse, whose semi-minor axis is equal to the width of the plume
    (i.e. the radius of the cylinder) and whose semi-major axis is determined by the
    viewing geometry.
    
    * r - radius of the plume (of the cylinder)
    * theta0 - scan angle of the centre of the plume (in degrees)
    * h_dist - the horizontal distance of the scanner from the centre of the plume
    * v_angle - the angle (in degrees) from horizontal of the scan (e.g. if you were 
                level with the plume then this will be zero).
    * theta - numpy array of scan angles.
    
    Returns: array of path lengths through the plume (one for each scan angle). Scan angles
             that did not intersect the plume will have a path length of zero.
    
    ---------------------------------------------------------------------------------------
    A few doctests...
    
    when v_angle=0 then we are dealing with a circle
    
    >>> p = length_of_scan_path(1.0, 90.0, 10.0, 0.0, numpy.array([90.0]))
    >>> print p[0] == 2.0
    True
    >>> p = length_of_scan_path(1.0, 90.0, 10.0, 0.0, numpy.array([0.0]))
    >>> print p[0] == 0.0
    True
    >>>
    
    when v_angle=60 then we expect the semi-major axis to be 2*r
    
    >>> p = length_of_scan_path(1.0, 90.0, 10.0, 60.0, numpy.array([90.0]))
    >>> print round(p[0],4) == 4.0
    True
    >>> p = length_of_scan_path(1.0, 90.0, 10.0, 60.0, numpy.array([0.0]))
    >>> print p[0] == 0.0
    True
    >>>
    
    #when theta is such that the line is a tangent, we expect to get zero back
    #>>> p = length_of_scan_path(1.0, 90.0, 10.0, 0.0, numpy.asin
    
    ---------------------------------------------------------------------------------------
    """
    #ellipse: x^2/a^2 + y^2/b^2 = 1
    #ellipse params are constrained by the viewing angle
    a = r    
    b = r / math.cos(math.radians(v_angle))
    
    #transform h_dist into the plane of the ellipse
    h_dist = h_dist / math.cos(math.radians(v_angle))
    
    #each scan line can be defined by two points on the line - here we use the 
    #x and y axis intercepts as the points
    x1 = numpy.zeros_like(theta)
    y1 = numpy.ones_like(theta) * (-h_dist)
    
    x2 = h_dist * numpy.tan(numpy.radians(theta - theta0))
    y2 = x1.copy()
    
    #use the equation of a straight line x = cy + d we do it 
    #this way around because it is possible that theta==theta0, which 
    #leads to x1==x2 and an infinite gradient. However, h_dist is never
    #0 and so dx/dy will never be infinite.
    c = (x2 - x1) / (y2 - y1)
    d = x1 - (c * y1)
    
    #substitute the expression for x from the line into the equation of the ellipse
    #and solve using the quadratic equation.
    #quadratic equation -> x = (-B +- sqrt(B^2 - 4AC))/(2A)
    A = (a ** 2 / b ** 2) + c ** 2
    B = 2.0 * c * d
    C = d ** 2 - a ** 2

    #mask out scan angles that did not intersect the plume (they have imaginary roots)        
    discriminant = B ** 2 - 4 * A * C
    invalid_values = discriminant <= 0
    
    A = numpy.ma.masked_where(invalid_values, A, copy=False)
    B = numpy.ma.masked_where(invalid_values, B, copy=False)
    C = numpy.ma.masked_where(invalid_values, C, copy=False)
    c = numpy.ma.masked_where(invalid_values, c, copy=False)
    d = numpy.ma.masked_where(invalid_values, d, copy=False)

    #solve for the intersection points of the lines with the ellipse    
    cy1 = (-B + numpy.sqrt(B ** 2 - 4.0 * A * C)) / (2.0 * A)
    cy2 = (-B - numpy.sqrt(B ** 2 - 4.0 * A * C)) / (2.0 * A)
        
    cx1 = c * cy1 + d
    cx2 = c * cy2 + d
    
    #return the distances between the intersection points (i.e. the path lengths
    #and fill in the values that did not intersect with zeros.
    return numpy.ma.filled(numpy.ma.sqrt((cx2 - cx1) ** 2 + (cy2 - cy1) ** 2), 0.0)
    


def estimate_plume_radius(angles, ydata, h_dist, v_angle, subplot=None):
    """
    Assuming a cylindrical plume, this function finds the plume radius and 
    position (scan angle to the centre) that best fits the supplied data.
    It works by doing a least squares fit of linearly scaled path lengths
    through the plume for different radii to the data. In this way it is 
    assuming that the plume is uniform and that changes in column amounts
    are only due to different path lengths.
    
    * angles - numpy array of scan angles
    * ydata - numpy array of column amounts or ratios or similar
    * h_dist - the horizontal distance of the scanner from the centre of the plume
    * v_angle - the angle (in degrees) from horizontal of the scan (e.g. if you were 
                level with the plume then this will be zero).
    * subplot - a matplotlib Axes object that you want the best fit function to be
                plotted into, or None if you don't want anything plotted.
                
    Returns: a tuple of best fit parameters (radius, theta0) where theta0 is the scan
             angle towards the centre of the plume.
    
    ---------------------------------------------------------------------------------------
    A few doctests...
    
    >>> ang = numpy.array(range(45, 136))
    >>> so2 = length_of_scan_path(200.0, 90.0, 1000.0, 45, ang)
    >>> r,theta0 = estimate_plume_radius(ang, so2, 1000.0, 45)
    >>> print round(r,3) == 200.0
    True
    >>> print round(theta0,3) == 90.0
    True
    
    """

    # guess some fit parameters - we use the width at half height as an estimate of radius
    ymax = ydata.max()

    #print "weighted av = ", numpy.average(angles,weights=numpy.abs(ydata))

    angle_at_half_height = angles[numpy.abs(ydata - (ymax / 2.0)).argmin()]
       
    theta0_guess = numpy.average(angles,weights=numpy.abs(ydata))
    
    variance = numpy.dot(numpy.abs(ydata), (angles-theta0_guess)**2)/numpy.abs(ydata).sum()  # Fast and numerically precise
    
    r_guess_angle = math.sqrt(variance)
    r_guess = h_dist * math.tan(math.radians(r_guess_angle))
    
    #print "theta0 guess = ", theta0_guess
    
    #r_guess = h_dist * math.tan(math.radians(math.fabs(theta0_guess - angle_at_half_height)))
    #print "r_guess = ",r_guess
    scaling_guess = ymax / (2.0 * (r_guess / math.cos(math.radians(v_angle))))

    #put guess params into an array ready for fitting
    p0 = numpy.array([r_guess, theta0_guess, scaling_guess])

    #define the path length function and associated error function
    fitfunc = lambda p, x: p[2] * length_of_scan_path(p[0], p[1], h_dist, v_angle, x)
    errfunc = lambda p, x, y: fitfunc(p, x) - y
   
    # do the fitting
    p1, success = scipy.optimize.leastsq(errfunc, p0, args=(angles, ydata))

    if success not in (1, 2, 3, 4):
        raise RuntimeError, "Could not fit cylindrical plume to data."

    if subplot is not None:
        subplot.plot(angles, ydata, 'g+', label='data')
        
        #plot the ellipse fit at higher angle resolution
        high_res_ang = numpy.linspace(numpy.min(angles), numpy.max(angles), 1000)
        subplot.plot(high_res_ang, fitfunc(p1, high_res_ang) , 'r-', label='fit')
    
    return (p1[0], p1[1])



def find_residuals(x, y, x_err, y_err, func):
    """
    Calculates the absolute distances of a data series from a trend line.
    However, any distances that are within the error bars are set 
    to zero. This makes it easy to see if a trend line fits some data within
    experimental error. 
    
    * x - numpy array of x data
    * y - numpy array of corresponding y data
    * x_err - absolute error in x direction (x +- x_err)
    * y_err - absolute error in y direction (y +- y_err)
    * func - a Python function that describes the trendline such that y_trend = func(x)
             where x is a single x value.
    """
    #convert func to numpy ufunc
    fit_func = numpy.frompyfunc(func, 1, 1)
    
    x_mins = x - x_err
    x_maxs = x + x_err
    
    #evaluate the distance between the data and the trendline at the extremes of
    #the x error.
    resid_mins = numpy.abs(y - fit_func(x_mins))
    resid_maxs = numpy.abs(y - fit_func(x_maxs))
    
    #get the smallest distance of each point from the fit line
    resid = numpy.minimum(resid_mins, resid_maxs)
    
    #set residuals that are within the error bars to zero
    resid[numpy.where(resid < y_err)] = 0
    
    return resid



def fullscans(filename, conv=1.0):
    """
    Generator that take the gradient of the angle column, creates 2 boolean array of true and 
    false based on positive or negative gradient offset the 2 array by one then
    assign inflection points where the arrays are not equal. The generator yield
    a dictionary "d" that contain each subset corresponding to a fullscan. 
    """  
    times = []
    angle = []
    angle_err = []
    col_amount = []
    col_amount_err = []
    int_times = []
    n_coadds = []
    
    with open(filename,"r") as ifp:
        #load values and errors from file
        for line in ifp:
            if line.startswith("#") or line.isspace() or line == "":
                continue            
            words = line.split()
            try:
                times.append(datetime.datetime.strptime(words[0].lstrip().rstrip(),"%d/%m/%Y_%H:%M:%S.%f"))
            except ValueError:
                times.append(datetime.datetime.strptime(words[0].lstrip().rstrip(),"%d/%m/%Y_%H:%M:%S"))
            angle.append(float(words[2].lstrip().rstrip()))
            angle_err.append(float(words[3].lstrip().rstrip()))
            col_amount.append(conv*float(words[4].lstrip().rstrip()))
            col_amount_err.append(conv*float(words[5].lstrip().rstrip()))
            int_times.append(float(words[6].lstrip().rstrip()))
            n_coadds.append(int(words[7].lstrip().rstrip()))
            
    
    #convert to numpy arrays
    times = numpy.array(times)
    angle = numpy.array(angle)
    angle_err = numpy.array(angle_err)
    col_amount = numpy.array(col_amount)
    col_amount_err = numpy.array(col_amount_err)
    int_times = numpy.array(int_times)
    n_coadds = numpy.array(n_coadds)
    
    for i in processing.split_by_scan(angle,times, angle_err,col_amount,col_amount_err, int_times, n_coadds):
        d = {}
        d["angle"] = i[0]
        d["times"] = i[1]
        d["angle_error"] = i[2]
        d["so2_col_amount"] = i[3]
        d["so2_col_amount_error"] = i[4]
        d["int_times"] = i[5]
        d["coadds"] = i[6]
        
        yield d
 
    


def calculate_bkgd(angles, col_amount, angle_err, col_err, fit_order=2, subplot_limits=None, subplot_calc=None):
    """
    Calculates a polynomial fit to the background for a single scan. This is done
    by fitting a Gaussian to the scan, and treating all points 2.5*sigma away from the
    mean as background points. Points which cannot be fitted within error are then removed
    one by one and the fitting is repeated each time.
    
    * angles - numpy array of scan angles
    * col_amount - numpy array of column amounts or ratios or similar
    * angle_err - numpy array of absolute errors in scan angles
    * col_err - numpy array of absolute errors in the column amounts
    * fit_order - order of polynomial that you want to fit
    * subplot - a matplotlib Axes object that you want the background fit function to be
                plotted into, or None if you don't want anything plotted.
    
    Returns: a numpy array of background levels (with the same length as the supplied 
             col_amount array).
    """
    
    #Do the Gaussian fitting
    gfit_params = processing.fit_gaussian(angles, col_amount)
       
    #mask the angles and col_amounts array to only expose the out of plume sections
    a = 2.5
    is_lower_bkgd = angles < (gfit_params.mean - a * gfit_params.sigma)
    is_higher_bkgd = angles > (gfit_params.mean + a * gfit_params.sigma)
    bkgd_mask = numpy.logical_or(is_lower_bkgd, is_higher_bkgd)
    bkgd_idxs = numpy.where(bkgd_mask)
    
    #fit a polynomial through the background
    ignored_pts = []
    
    while True:
        bkgd_angles = angles[bkgd_idxs]
        bkgd_col_amounts = col_amount[bkgd_idxs]
        bkgd_angle_err = angle_err[bkgd_idxs]
        bkgd_col_err = col_err[bkgd_idxs]
        
        polyfit_params = numpy.polyfit(bkgd_angles, bkgd_col_amounts, fit_order)
        bkgd_function = numpy.poly1d(polyfit_params)
             
        #calculate the residuals - ignoring those that are within error
        resid = find_residuals(bkgd_angles, bkgd_col_amounts, bkgd_angle_err, bkgd_col_err, bkgd_function)
        
        if len(numpy.nonzero(resid)[0]) == 0:
            break
        
        #remove the point with the greatest residual from the fit
        idx_to_skip = bkgd_idxs[0][numpy.argmax(resid)]
        ignored_pts.append((angles[idx_to_skip], col_amount[idx_to_skip], angle_err[idx_to_skip], col_err[idx_to_skip]))
        bkgd_idxs = (bkgd_idxs[0][numpy.where(bkgd_idxs[0] != idx_to_skip)],)
    
    print "Ignored", len(ignored_pts), "points from background fitting."
   
    #if we have been passed a subplot object - then plot the fit into it     
    if subplot_calc is not None:
        subplot_limits.errorbar(angles, col_amount, col_err, angle_err, '+')
        globals()["d"] = DraggableLine(gfit_params.mean - a * gfit_params.sigma, subplot_limits)
        globals()["d2"] = DraggableLine(gfit_params.mean + a * gfit_params.sigma, subplot_limits)
        #subplot_limits.plot([gfit_params.mean - a * gfit_params.sigma]*2,subplot_limits.get_ylim(),'b-', linewidth=2)
        #subplot_limits.plot([gfit_params.mean + a * gfit_params.sigma]*2,subplot_limits.get_ylim(),'g-', linewidth=2)
        subplot_limits.grid()
         
        not_bkgd_idxs = numpy.where(numpy.logical_not(bkgd_mask))
        not_bkgd_angles = angles[not_bkgd_idxs]
        not_bkgd_col_amounts = col_amount[not_bkgd_idxs]
        not_bkgd_col_err = col_err[not_bkgd_idxs]
        not_bkgd_angle_err = angle_err[not_bkgd_idxs]
        
        ignored = zip(*ignored_pts)
        if ignored:
            subplot_calc.errorbar(ignored[0], ignored[1], ignored[3], ignored[2], 'r+')
        
        #subplot_calc.errorbar(not_bkgd_angles,not_bkgd_col_amounts, not_bkgd_col_err, not_bkgd_angle_err,'g+')        
        subplot_calc.errorbar(bkgd_angles, bkgd_col_amounts, bkgd_col_err, bkgd_angle_err, 'b+')
        subplot_calc.plot(angles, [bkgd_function(i) for i in angles], 'm-', linewidth=2)
    
    #return the calculated background for all angles
    return numpy.array([bkgd_function(i) for i in angles])



def integrate(angles, data, radius, h_dist, v_angle, theta0):
    """
    
    
    ---------------------------------------------------------------------------------------
    A few doctests...
    
    >>> ang = numpy.array(range(45, 136))
    >>> so2 = length_of_scan_path(200.0, 90.0, 1000.0, 45, ang)
    >>> i = integrate(ang, so2, 200.0, 1000.0, 45.0, 90.0)
    >>> print round(i, 4) == round(200**2 * math.pi, 4)
    True
  
    """
    #TODO - data needs to be converted to something per m^2 rather than whatever it is in currently
    #TODO - error analysis!
    
    #the integration assumes the angles are in increasing order so sort them here
    c = numpy.rec.fromarrays([angles, data])
    c.sort() #sort based on angles
    angles = c.f0
    data = c.f1

        
    #calculate the path length through the plume for each of the measurements
    measured_dists = cPlume.get_path_length(radius, theta0, h_dist, v_angle, angles)
    
    #crop data to only be working inside the plume
    pts_inside_plume = numpy.where(measured_dists > 0)
    angles = angles[pts_inside_plume]
    data = data[pts_inside_plume]
    measured_dists = measured_dists[pts_inside_plume]
    
    #if there are no points inside the plume - return 0
    if len(data) == 0:
        return 0.0
    
    
    #calculate the areas of the circle sections between each scan line
    #for chord area calculation see: http://mathworld.wolfram.com/Chord.html
    r = h_dist * numpy.sin(numpy.radians(angles - theta0))    
    chord_areas = numpy.zeros((len(angles) + 1,), dtype="float")
    chord_areas[-1] = (math.pi * radius ** 2)
    chord_areas[:-1] = chord_areas[-1] - ((radius ** 2 * numpy.arccos(r / radius)) - (r * numpy.sqrt(radius ** 2 - r ** 2)))
    
    #scan areas are equal to the difference between chord areas
    scan_areas = numpy.zeros_like(chord_areas)
    scan_areas[0] = chord_areas[0]
    scan_areas[1:] = chord_areas[1:] - chord_areas[:-1]
       
    #calculate the averaged concentraions
    concs = data / measured_dists
    
    av_concs = numpy.zeros_like(scan_areas)
    av_concs[0] = concs[0] #don't bother averaging at the edges
    av_concs[-1] = concs[-1] #don't bother averaging at the edges
    av_concs[1:-1] = (concs[:-1] + concs[1:]) / 2.0 
    
    
    #calculate the integral
    return numpy.sum(av_concs * scan_areas)
 


def generate_random_angles_array(times, angles, angle_err, n_iter, n_coadds, int_times ):
    
    #convert times to seconds
    times = numpy.array([date2secs(t) for t in times])
    
    #calculate the minimum distance that points can be apart - integration time * number of coadds
    min_spacings = numpy.zeros_like(angles)
    min_spacings[1:] = ((angles[1:] - angles[:-1]) / (times[1:] - times[:-1])) * n_coadds[1:] * int_times[1:]
    
    #generate an array of uniformly distributed random numbers between 0 and 1
    rand_nums = numpy.random.random(size=(n_iter,len(angle_err)))
    
    rand_angle_errs = numpy.zeros_like(rand_nums)
    
    #first column can take full range of error between +/- angle_err
    rand_angle_errs[:,0] = ((2*angle_err[0]) * rand_nums[:,0]) - angle_err[0]
    
    for i in range(1,len(angle_err)):
        a1 = angles[i] - (angles[i-1] + rand_angle_errs[:,i-1]) - min_spacings[i] #should be min_spacings[i] (not i-1) because need int_time of current point
        
        #if the minus error is less than the negative of the plus error then there is a problem - 
        #the point can't be fitted in with the required gaps
        if numpy.any(-a1 > angle_err[i]) :
            print a1[numpy.where(-a1 > angle_err[i])]
            print angle_err[i]
            
            if angles[i] == angles[i-1]:
                print "equal angles is the problem!"
            raise ValueError, "No space left for data point!"
        
        a = numpy.minimum(a1,angle_err[i])
        
        rand_angle_errs[:,i] = ((angle_err[i]+a)*rand_nums[:,i]) - a
    
    return rand_angle_errs

#This is the old version of this function which integrates the entire scan
#(not just the bit between the plume-limits). This causes problems if the scan
#hits the ground or has anomolous points in the background, since these will
#be included in the flux estimate.
#def calculate_ica_trapz(scan):
#    
#    #subtract the background from the col_amounts
#    col_amounts = scan.col_amount - scan.bkgd_func(scan.angles)
#    
#    #reproject the column amounts given the inclination of the scan
#    corrected_col_amounts = math.cos(math.radians(scan.v_angle)) * col_amounts
#    
#    #calculate distance between measurements assuming dx = r * theta
#    d_theta = numpy.radians(numpy.abs(scan.angles[1:]-scan.angles[:-1]))
#    dx = scan.h_dist * d_theta
#    
#    #calculate the error in the distances
#    #first the random error in the angle
#    d_theta_err = numpy.sqrt(numpy.radians(scan.angle_err[1:])**2 + numpy.radians(scan.angle_err[:-1])**2)
#    
#    #systematic error in the distance
#    dist_error = abs(0.5 * scan.h_dist_err * math.cos(math.radians(scan.v_angle)) * numpy.sum(d_theta * (col_amounts[:-1] + col_amounts[1:])))
#    
#    #systematic error in v_angle
#    v_angle_error = abs(-0.5 * math.radians(scan.v_angle_err) * math.sin(math.radians(scan.v_angle)) * numpy.sum(d_theta * (col_amounts[:-1] + col_amounts[1:])))
#    
#    #random errors in angle and colulm amount
#    col_amt = dx * ((corrected_col_amounts[:-1] + corrected_col_amounts[1:])/2.0)
#
#    delta_col_amt = col_amt * numpy.sqrt((d_theta_err/d_theta)**2 + (numpy.sqrt(scan.col_err[:-1]**2 + scan.col_err[1:]**2) / (corrected_col_amounts[:-1] + corrected_col_amounts[1:]))**2)
#    
#    #anywhere that d_theta was zero - set the error to zero (otherwise end up with NaNs in array
#    delta_col_amt[numpy.where(d_theta==0)]=0
#    
#    random_errors = math.sqrt(numpy.sum(delta_col_amt**2))
#    
#
#    
#    #use the trapezoid rule to integrate
#    print "setting scan.ica"
#    scan.ica = numpy.sum(col_amt)
#    
#    scan.ica_err = dist_error + v_angle_error + random_errors
#        
#    min_time = min(scan.times)
#    max_time = max(scan.times)
#    scan.ica_time = min_time + ((min_time - max_time) / 2)
#    scan.ica_time_err = datetime.timedelta(seconds=0)
    
    
    
def calculate_ica_trapz(scan):
    
    times, angles, col_amounts, angle_err, col_err, coadds, int_times = scan.get_bkgd_subtracted_plume_points()
    
    #if there is no plume region selected then set everything to 0
    if len(times) == 0:
        scan.ica = 0.0
        scan.ica_err = 0.0
        min_time = min(scan.times)
        max_time = max(scan.times)
        scan.ica_time = min_time + ((min_time - max_time) / 2)
        scan.ica_time_err = 0.0
    
    #reproject the column amounts given the inclination of the scan
    corrected_col_amounts = math.cos(math.radians(scan.v_angle)) * col_amounts
    
    #calculate distance between measurements assuming dx = r * theta
    d_theta = numpy.radians(numpy.abs(angles[1:]-angles[:-1]))
    dx = scan.h_dist * d_theta
    
    #calculate the error in the distances
    #first the random error in the angle
    d_theta_err = numpy.sqrt(numpy.radians(angle_err[1:])**2 + numpy.radians(angle_err[:-1])**2)
    
    #systematic error in the distance
    dist_error = abs(0.5 * scan.h_dist_err * math.cos(math.radians(scan.v_angle)) * numpy.sum(d_theta * (col_amounts[:-1] + col_amounts[1:])))
    
    #systematic error in v_angle
    v_angle_error = abs(-0.5 * math.radians(scan.v_angle_err) * math.sin(math.radians(scan.v_angle)) * numpy.sum(d_theta * (col_amounts[:-1] + col_amounts[1:])))
    
    #random errors in angle and colulm amount
    col_amt = dx * ((corrected_col_amounts[:-1] + corrected_col_amounts[1:])/2.0)

    delta_col_amt = col_amt * numpy.sqrt((d_theta_err/d_theta)**2 + (numpy.sqrt(col_err[:-1]**2 + col_err[1:]**2) / (corrected_col_amounts[:-1] + corrected_col_amounts[1:]))**2)
    
    #anywhere that d_theta was zero - set the error to zero (otherwise end up with NaNs in array
    delta_col_amt[numpy.where(d_theta==0)]=0
    
    random_errors = math.sqrt(numpy.sum(delta_col_amt**2))
    

    
    #use the trapezoid rule to integrate
    print "setting scan.ica"
    scan.ica = numpy.sum(col_amt)
    
    scan.ica_err = dist_error + v_angle_error + random_errors
        
    min_time = min(scan.times)
    max_time = max(scan.times)
    scan.ica_time = min_time + ((min_time - max_time) / 2)
    scan.ica_time_err = datetime.timedelta(seconds=0)
    
     
def calculate_ica_cyl(scan, n_iter):
    
    times, raw_angles, col_amount, angle_err, col_err, coadds, int_times = scan.get_bkgd_subtracted_plume_points()
    
    if len(raw_angles) < 5:
        #if there are less than 5 points in the plume, then the background is probably
        #badly defined - don't bother trying to calculate the ICA, just set it to zero
        scan.ica = 0
        scan.ica_err = 0
        
        min_time = min(scan.times)
        max_time = max(scan.times)
        scan.ica_time = min_time + ((min_time - max_time) / 2)
        scan.ica_time_err = abs(scan.ica_time - min_time)
        return 
    
    if n_iter > 1:
        #create arrays of random errors (normally distributed with std_dev=1.0)
        #rnd_nums_angles = numpy.random.normal(size=(n_iter,len(raw_angles)))
        #rnd_nums_angles = rnd_nums_angles*angle_err
        rnd_nums_angles = generate_random_angles_array(times, raw_angles, angle_err, n_iter, coadds, int_times)
        
        rnd_nums_col_amounts = numpy.random.normal(size=(n_iter,len(col_amount)))
        rnd_nums_col_amounts = rnd_nums_col_amounts*col_err
        
        #add the random errors to our data
        angles = rnd_nums_angles + raw_angles
        col_amount = rnd_nums_col_amounts + col_amount
        
        #it is unphysical for two of the values to be swapped in angle due to error on the angle i.e. the 
        #order that the data is in must be preserved. Any set of angles which contain two or more points whose
        #order has been swapped due to the addition of random error are replaced here by a new set with randomly
        #generated errors
        for i in range(angles.shape[0]):
            j=0
            while numpy.any(angles[i] != numpy.sort(angles[i])): 
                print "Row ",i," contains an overlap - regenerating random errors."
                new_rnd_nums = numpy.random.normal(size=(1,len(angles[i])))
                new_rnd_nums = new_rnd_nums*(2*angle_err)-angle_err
                angles[i] = new_rnd_nums + raw_angles[i]
                j+=1
                if j>1000:
                    print "Too much chance of value overlap. I give up!"
                    return
            
    else:
        #the arrays need to be 2D even if they only have one row
        angles = raw_angles.reshape((1,len(raw_angles)))  
        col_amount = col_amount.reshape((1,len(col_amount)))

        
    #make a first guess at fit parameters
    ymax = col_amount.max(axis=1)

    theta0_guess = numpy.average(angles, axis=1, weights=numpy.abs(col_amount))
    
    #TODO - do variance calculation as an array operation
    variance = numpy.zeros((n_iter,1), dtype=float)
    for i in range(n_iter):                         
        variance[i] = numpy.dot(numpy.abs(col_amount[i]), (angles[i]-theta0_guess[i])**2)/numpy.abs(col_amount[i]).sum()  # Fast and numerically precise
            
    r_guess_angle = numpy.sqrt(variance)
    r_guess = (scan.h_dist * numpy.tan(numpy.radians(r_guess_angle))).flatten()


    scaling_guess = ymax / (2.0 * (r_guess / math.cos(math.radians(scan.v_angle))))

    #put guess params into an array ready for fitting
    p0 = numpy.array([r_guess, theta0_guess, scaling_guess]).swapaxes(1,0)
    
    #define the path length function and associated error function
    fitfunc = lambda p, x: p[2] * cPlume.get_path_length(p[0], p[1], scan.h_dist, scan.v_angle, x)
    errfunc = lambda p, x, y: fitfunc(p, x) - y     
    
    # do the fitting
    icas = numpy.zeros((n_iter,1), dtype='float')
    p1 = numpy.zeros_like(p0)
    measurement_times = []
    max_fev = 200*len(angles[0])+200 #increase the max number of iterations for the fitting (otherwise sometimes it doesn't find the fit)
    for i in range(n_iter):
        ang = angles[i]
        so2 = col_amount[i]
        p1[i], success = scipy.optimize.leastsq(errfunc, p0[i], args=(ang, so2), maxfev=max_fev)
        
        if success not in (1, 2, 3, 4):
            #for scans that don't have many points inside the plume, this quite often happens
            #the shape of the plume is not well enough constrained to fit the ellipse function
            #to it (at least not within the allowed number of iterations. I can't think up a 
            #smarter thing to do other than to give up on the error analysis and recalculate the
            #point with no errors and warn the user that this has happened.
            print "You can ignore the above warning.....I've got this."
            if n_iter == 1:
                #uh oh! we are in danger of infinite recursion here. Set the ica to zero
                #and return
                scan.ica = 0
                scan.ica_err = 0
                scan.ica_time_err = datetime.timedelta(seconds=0)
            else:
                calculate_ica(scan,1)
            return
            

        icas[i] = integrate(ang, so2, p1[i][0], scan.h_dist, scan.v_angle, p1[i][1])
        
        idx_theta0 = numpy.abs(ang - p1[i][1]).argmin()
        measurement_times.append(scan.times[idx_theta0])
       
    #update the scan parameters
    scan.ica = numpy.mean(icas)
    scan.ica_err = numpy.std(icas) 
    scan.error_calc_iters = n_iter

    av_theta0 = p1[:,1].mean()
    av_radius = p1[:,0].mean()
    
    av_radius_degrees = math.degrees(math.atan(av_radius/scan.h_dist))

    scan.ica_time = scan.angle2time(av_theta0)

    if n_iter >1:
        scan.ica_time_err = abs(scan.ica_time - scan.angle2time(av_theta0 - av_radius_degrees))
    else:
        scan.ica_time_err = datetime.timedelta(seconds=0)
